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ABSTRACT 

We  apply  the  boundary  element  method  (BEM)  to  study  seismic  wave  radiation  from  explosions 
in  cylindrical  cavities.  This  approach  simulates  wave  propagation  by  replacing  the  boundary 
between  the  cavity  wall  and  the  surrounding  formation  by  a  distribution  of  fictitious  sources 
that  produce  the  same  composite  wavefield  as  the  original  model.  The  cavity  wall  is  discretized 
into  elements  in  the  numerical  implementation,  and  the  fictitious  source  amplitudes  on  each 
element  are  determined  by  solving  boundary  condition  equations.  In  order  to  accomplish  this, 
the  interactions  of  all  elements  must  be  computed  as  well  as  the  wavefields  radiating  from 
the  sources  to  the  receiver  positions.  A  parallel  computer  implementation  of  the  algorithm 
is  therefore  advantageous,  since  the  frequency  domain  computations  are  easily  parallelized  for 
rapid  solution  of  the  problem.  We  also  compute  discrete  wavenumber  synthetic  seismograms 
at  regional  distances  from  the  cavity  source,  using  an  equivalent  moment  tensor  model  of  the 
explosion  in  the  cylindrical  cavity. 

Results  for  the  radiation  patterns  show  that  at  large  aspect  ratios  (the  ratio  of  cavity  length 
to  cavity  diameter),  the  source  has  a  strong  variation  in  signal  strength  at  high  frequency. 
However,  at  lower  frequencies  typical  of  regional  wave  propagation,  the  P-wave  radiation  pattern 
is  more  nearly  isotropic,  and  the  S-wave  amplitude  is  not  so  large.  In  addition,  we  find  the 
interesting  result  that  the  radiation  pattern  of  a  large  aspect  ratio  cavity  in  a  transversely 
isotropic  medium  can  be  more  like  a  classical  point  source  than  an  a  similar  isotropic  example. 
Incorporation  of  the  cavity  source  into  simulations  of  regional  wave  propagation  using  discrete 
wavenumber  techniques  shows  that  even  though  the  shear  wave  is  relatively  weak  at  the  source, 
it  can  still  lead  to  the  generation  of  a  significant  transverse  component  Lg  wave  at  regional 
distances. 
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OBJECTIVE 


The  primary  goal  of  the  work  reported  in  this  paper  is  to  examine  the  influence  of  source 
cavity  geometry  on  radiation  of  seismic  waves  at  regional  distances.  Specifically,  by  examining 
radiation  patterns  of  these  sources,  we  seek  to  determine  whether  or  not  these  sources  can 
explain  variations  in  S  and  P-wave  energy  observed  in  field  data,  and  whether  or  not  they  can 
generate  significant  SH-waves.  We  therefore  simulate  the  radiation  patterns  of  explosion  sources 
located  in  cylindrical  cavities  that  model  tunnels  using  the  boundary  element  method  (BEM).  In 
addition,  we  examine  the  behavior  of  such  sources  in  representative  models  of  the  earth’s  crust 
using  the  discrete  wavenumber  method  and  have  computed  results  for  these  sources  in  3-D  earth 
models  using  finite  difference  methods.  Results  from  these  numerical  experiments  show  that 
the  radiation  patterns  from  explosion  sources  located  in  tunnels  may  not  produce  significantly 
strong  shear  waves.  If  the  rock  formations  surrounding  the  borehole  are  transversely  isotropic 
(axis  of  symmetry  parallel  to  the  cavity  axis),  the  BEM  results  show  that  the  quasi-shear  wave 
produced  in  this  case  can  actually  be  weaker  than  in  the  isotropic  case.  Discrete  wavenumber 
seismograms  show  that  even  though  the  shear  wave  is  somewhat  weaker  than  the  compressional 
wave  at  the  source,  it  still  produces  a  strong  Lg  wave  at  regional  distances  with  some  amplitude 
variations  as  a  function  of  observation  azimuth. 

RESEARCH  ACCOMPLISHED 


Introduction 

The  radiation  patterns  of  seismic  sources  are  very  important  in  utihzing  seismic  signals  to 
distinguish  between  nuclear  explosions  and  natural  earthquake  events.  The  standard  model 
of  the  nuclear  explosion  source,  the  spherically  symmetric  pressure  field  corresponding  to  a 
moment  tensor  model  with  three  equal,  mutually  perpendicular  dipoles,  generates  only  an 
isotropic  compressional  wave  disturbance.  However,  any  deviation  from  this  ideal  configuration 
can  significantly  change  the  radiation  pattern  of  the  source.  For  example,  Zhao  and  Harkrider 
(1992)  show  that  when  the  source  is  located  off-center  in  a  spherical  cavity  it  will  generate 
significant  shear  wave  energy.  Their  work  considers  a  source  located  in  a  spherical  solid,  which 
is  in  turn  embedded  within  another  infinite  sohd,  representing  a  fully  tamped  explosion. 

Another  way  of  altering  the  radiation  of  the  source  is  to  change  the  shape  of  the  cavity. 
Glenn  et  al.  (1985,  1986),  Rial  and  Moran  (1986)  and  Glenn  and  Rial  (1987)  considered  the 
radiation  patterns  of  sources  located  in  elhpsoidal  and  cylindrical  cavities.  Using  a  combination 
of  analytical  and  numerical  methods,  they  obtained  results  suggesting  that  the  P-wave  radiation 
pattern  would  be  strongest  perpendicular  to  the  cavity  axis  and  that  a  strong  shear  wave  would 
be  generated  as  cavity  aspect  ratio  (ratio  of  length  to  diameter)  became  large.  However,  Stevens 
et  al.  (1991)  performed  non-linear  numerical  simulations  that  suggested  that  at  lower  frequencies 
more  typical  of  regional  wave  propagation,  the  compressional  wave  radiation  pattern  from  a 
cylindrical  cavity  would  be  relatively  isotropic  and  only  a  small  S-wave  would  be  generated. 
More  recently,  Ben-Menahem  and  Mikhailov  (1995)  derived  approximate  analytic  solutions 
for  radiation  from  cylindrical  cavities.  In  contrast  to  the  results  of  Stevens  et  al.  (1991), 
they  predict  a  comparatively  large  S-wave  amphtude  for  low  frequencies.  Hence,  there  are 
significantly  different  results  for  cylindrical  cavity  radiation  patterns  available  in  the  literature. 

We  therefore  continue  this  line  of  investigation  by  applying  the  boundary  element  method 


375 


(BEM)  to  simulate  the  radiation  of  seismic  waves  from  a  cylindrical  cavity  representing  a  source 
located  in  a  vertical  shaft  or  a  horizontal  tunnel.  The  BEM  is  well  suited  to  problems  such 
as  this,  where  the  model  combines  small  features  (the  cavity)  with  large  features  (propagation 
distances  on  the  order  of  10  km  or  more)  (Bouchon,  1993;  Dong,  1993;  Dong  et  al.,  1995; 
Dong  and  Toksoz,  1995).  After  a  brief  description  of  this  method,  we  show  radiation  patterns 
obtained  with  the  BEM  that  are  similar  to  those  obtained  by  Stevens  et  al.  (1991).  They  show 
a  weak  variation  of  P-wave  amplitude  with  direction  and  a  comparatively  small  shear  wave 
at  frequencies  characteristic  of  regional  wave  propagation.  After  presenting  these  results  and 
demonstrating  the  effects  of  locating  the  cavity  in  a  transversely  isotropic  medium,  we  consider 
the  effects  of  such  a  source  in  a  plane  layered  model  of  the  earth.  The  source  is  applied  in 
a  discrete  wavenumber  code  by  developing  an  equivalent  moment  tensor  representation  of  the 
source.  Synthetic  seismograms  at  regional  distances  show  that  the  effect  of  the  cavity,  though 
subtle,  will  generated  SH-wave  motion  and  may  still  be  distinguishable  in  terms  of  an  azimuthal 
variation  in  amplitudes. 

Method 

The  BEM  has  recently  been  applied  to  synthesize  borehole  wave  propagation  in  experimental 
configurations  typical  of  exploration  geophysics  applications  (Bouchon,  1993;  Dong,  1993;  Dong 
et  al.,  1995;  Dong  and  Toksoz,  1995).  One  model  considered  in  this  earlier  work  is  a  semi¬ 
infinite  borehole,  where  the  seismic  source  is  apphed  at  or  near  the  end  of  the  borehole.  A 
straightforward,  logical  extension  of  this  approach  is  to  model  a  cylindrical  cavity  of  finite 
length,  such  as  the  source  cavities  that  might  be  utilized  in  testing  nuclear  devices.  Like  the 
exploration  simulations,  we  model  the  interior  of  the  cavity  as  an  acoustic  medium  and  also 
assume  symmetry  about  the  cavity  axis. 

The  essence  of  the  method  is  that  the  effects  of  the  diffracting  boundary  between  the  source 
cavity  and  the  surrounding  medium  can  be  reproduced  by  a  discrete  set  of  fictitious  secondary 
sources  (Bouchon,  1993).  Hence,  implementation  of  the  algorithm  begins  by  discretizing  the 
wall,  both  the  cylindrical  portion  and  the  flat  ends  capping  the  cavity.  The  magnitudes  of  the 
fictitious  sources  are  determined  by  enforcing  boundary  conditions  at  the  cavity  wall  (continuity 
of  radial  displacement  and  stress,  vanishing  of  shear  stress).  These  three  boundary  conditions 
are  sufficient  to  determine  the  sources,  since  for  each  element  there  is  a  volume  injection  source 
radiating  energy  into  the  acoustic  material  in  the  cavity,  and  a  vertical  and  radial  force  generat¬ 
ing  elastic  waves  in  the  rock  formation  outside  the  source  region.  The  determination  of  source 
magnitudes  is  accomplished  by  computing  the  displacement  field  generated  by  each  element 
incident  on  the  other  elements  using  a  discrete  wavenumber  method  to  evaluate  the  integrals 
over  horizontal  wavenumber.  Given  the  intense  computation  that  must  be  performed  in  order 
to  obtain  these  quantities,  it  is  advantageous  to  implement  the  algorithm  on  a  parallel  com¬ 
puter.  The  computations  for  each  frequency  are  independent,  so  we  assign  each  processor  on 
an  nCUBEl-2  computer  the  work  for  a  different  frequency. 

In  order  to  solve  the  resulting  matrix  of  boundary  condition  equations,  we  must  also  specify 
the  initial  field  generated  by  the  true,  primary  source.  A  straightforward  model  of  this  field, 
when  there  is  complete  decoupling  of  the  explosion,  is  to  assume  that  the  explosion  instanta¬ 
neously  pressurizes  the  cavity.  In  computing  our  frequency  domain  results,  we  therefore  model 
the  source  condition  as  a  step  function  of  constant  pressure  applied  to  the  cylindrical  cavity 
wall.  Earlier  versions  of  this  work  (Gibson  et  al.,  1994)  instead  used  a  volume  injection  point 
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Figure  1:  P-wave  amplitude  spectra  of  the  seismic  wavefields  radiated  by  a  source  in  a  cylindrical 
cavity  of  aspect  ratio  5.  Spectra  are  presented  for  three  different  propagation  directions  relative 
to  the  cavity  axis,  where  0°  is  parallel  to  the  axis.  The  directions  corresponding  to  3°  and  87° 
were  used  instead  of  0  and  90°  to  avoid  some  numerical  problems  associated  with  computations 
for  these  directions. 


source  located  on  the  axis  of  the  cavity.  The  differences  from  the  results  we  show  below  are 
caused  by  the  different  definition  of  the  initial  wavefields.  Although  for  relatively  small  cav¬ 
ity,  the  incident  wavefield  generated  by  the  point  source  will  strike  the  cavity  wall  essentially 
simultaneously,  the  amplitude  of  the  incident  field  will  decay  with  distance.  This  change  in 
amplitude  alters  radiation  patterns  somewhat  and  makes  the  patterns  sensitive  to  the  location 
of  the  source  within  the  cavity. 

Results 

Isotropic  medium 

We  applied  the  BEM  to  a  cavity  with  aspect  ratio  5  located  in  a  medium  with  a  P-wave  velocity 
of  4000  m/s,  S-wave  velocity  2200  m/s  and  density  2.2  gm/cm^.  These  parameters  will  show 
the  behavior  of  the  source  in  a  hard  rock  site.  Spectra  for  observation  directions  of  3°,  45°, 
and  87°  (0°  is  along  the  axis  of  the  cylindrical  cavity,  90°  is  perpendicular  to  it)  show  that  the 
corner  frequency  of  the  wavefields  depends  on  this  direction  (Figure  1).  Specifically,  the  corner 
frequency  is  higher  for  observation  directions  perpendicular  to  the  cavity.  These  results  are 
similar  to  those  of  Glenn  et  al.  (1985,  1986)  and  Stevens  et  al.  (1991). 

However,  most  of  the  amplitude  difference  in  these  wavefields  does  occur  at  high  frequencies. 
To  test  the  effects  that  might  be  seen  in  regional  seismograms,  we  computed  radiation  patterns 
at  lower  frequencies.  Radiation  patterns  of  both  P  and  S- waves  were  obtained  by  first  computing 
synthetic  seismograms  for  a  set  of  observation  points  at  a  constant  distance  of  10000  m  from 
the  source  and  convolving  the  frequency  domain  results  with  a  Ricker  wavelet  of  1  Hz. 

In  Figure  2,  we  show  the  resulting  radiation  patterns  for  cavities  80  m  long  with  varying 
aspect  ratio.  In  general,  the  changes  in  these  radiation  patterns  are  relatively  subtle,  but  some 
variations  are  evident.  As  aspect  ratio  decreases  from  5  to  values  near  1,  the  amplitude  of  the 
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ASPECT  RATIO  5  ASPECT  RATIO  2.5 


ASPECT  RATIO  1.25  ASPECT  RATIO  0.5 


Figure  2:  Radiation  patterns  for  explosion  sources  in  a  hard  rock  formation  (velocities  specified 
in  the  text).  Results  are  shown  for  several  values  of  the  cavity  aspect  ratio  to  demonstrate  the 
effect  of  variations  in  cavity  shape  on  far-field  seismic  wave  radiation. 


shear  wave  decreases  relative  to  the  compressional  wave,  and  the  compressional  wave  pattern 
becomes  more  isotropic.  When  the  aspect  ratio  is  less  than  one,  the  shear  wave  amplitude 
increases  once  again.  In  the  end  however,  the  changes  suggest  that  the  shear  wave  amplitude 
will  be  less  than  the  P-wave  at  frequencies  representative  of  regional  seismic  wave  propagation, 
a  conclusion  also  reached  by  Stevens  et  al.  (1991).  Ben-Menahem  and  Mikhailov  (1995),  on 
the  other  hand,  find  a  much  larger  S-wave  amplitude  using  an  approximate  analytic  solution 
that  maps  the  cylindrical  cavity  surface  onto  an  equivalent  spherical  cavity.  The  differences  in 
results  may  be  caused  by  this  mapping  for  a  cylindrical  cavity,  which  is  not  very  sphere-like. 

Transversely  isotropic  formation 

The  implementation  of  the  BEM  simulation  assumes  axial  symmetry  of  both  the  cavity  model 
and  the  wave  propagation.  A  transversely  isotropic  medium  can  be  incorporated  into  the 
modeling  as  long  as  the  axis  of  symmetry  of  the  formation  is  parallel  to  the  axis  of  symmetry 
of  the  cylindrical  cavity.  An  example  of  this  situation  would  be  a  medium  containing  aligned 
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ASPECT  RA  TIO  1.25  ASPECT  RA  TIO  0.5 


Figure  3:  Radiation  patterns  for  explosion  sources  in  a  fractured  hard  rock  formation.  This 
medium,  with  elastic  constants  shown  in  Table  1,  contains  cracks  aligned  perpendicular  to  the 
cavity  axis.  It  is  therefore  transversely  isotropic.  At  large  aspect  ratio,  the  S-wave  radiation  is 
weaker  than  in  the  isotropic  formation  (Figure  2). 


vertical  cracks,  a  material  that  is  effectively  anisotropic  with  a  horizontal  axis  of  symmetry  (e.g., 
Crampin,  1981).  If  the  axis  of  a  tunnel  is  horizontal  and  is  oriented  perpendicular  to  the  vertical 
cracks  in  such  a  transversely  isotropic  medium,  we  can  model  the  radiation  patterns  with  the 
current  implementation  of  the  BEM.  The  elastic  moduli  of  the  medium  can  be  computed  using 
the  theory  of  Hudson  (1980,  1981).  We  take  the  elastic  properties  of  the  uncracked  rock  to  be 
the  same  as  the  isotropic  formation  considered  above  and  set  the  crack  density  to  ^  =  na^  =  0.05 
(n= number  density  of  cracks,  a=radius  of  the  penny-shaped  cracks).  We  also  assume  that  the 
cracks  are  dry  (empty).  Table  1  displays  the  elastic  constants  of  both  the  isotropic  material 
and  the  composite  medium  with  ahgned  cracks. 

Radiation  patterns  for  the  cyhndrical  cavities  in  this  formation  are  displayed  in  Figure  3. 
The  results  are  somewhat  surprising  in  that  the  shear  waves  are  somewhat  weaker  than  in 
the  isotropic  medium  for  cavities  with  large  aspect  ratio.  Theoretical  results  show  that  a  point 
source  in  an  anisotropic  material  can  have  generate  a  very  large  shear  wave  (e.g.,  Ben-Menahem 
et  al.,  1991).  Since  the  ideal  point  source  explosion  in  an  isotropic  medium  will  generate  no 
shear  waves  at  all,  this  is  a  dramatic  change,  and  our  initial  intuition  might  suggest  that  shear 
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waves  should  also  be  large  for  the  cavity  sources.  However,  our  results  show  that  when  the 
explosion  is  located  in  the  cylindrical  cavity,  the  shape  of  the  cavity  alters  the  effective  moment 
tensor  corresponding  to  the  source  in  such  a  way  that  the  resulting  wave  field  resembles  the 
classic  isotropic  explosion  more  than  the  cylindrical  cavity  source  in  an  isotropic  medium! 

Regional  synthetic  seismograms 

Although  the  radiation  patterns  show  that  the  amplitude  of  the  shear  wave  is  not  too  large 
(Figures  2,  3),  the  effects  of  these  sources  at  regional  distances  is  not  necessarily  obvious.  In 
order  to  more  directly  predict  the  effects  of  such  source  within  layered  models  of  the  crust, 
we  computed  discrete  wavenumber  synthetic  seismograms  (e.g.,  Bouchon,  1982).  The  layer 
thickness  and  velocities  are  shown  in  Table  2.  The  receivers  were  located  on  the  free  surface  at 
a  distance  of  500  km  from  the  explosion  source,  and  the  full  range  of  azimuths  was  considered  to 
see  how  displacement  fields  would  depend  on  propagation  direction  relative  to  the  cavity  axis. 
The  explosion  source  itself  was  incorporated  within  the  modeling  by  estimating  the  moment 
tensor  producing  the  radiation  pattern  for  the  aspect  ratio  5  cavity  in  Figure  2  (tests  show 
that  the  moment  tensor  representation  does  not  depend  too  much  on  the  velocity  outside  the 
cavity). 

The  seismograms  show  that  the  vertical  component  seismograms  do  not  change  much  with 
azimuth,  and  they  look  essentially  the  same  as  the  seismograms  computed  for  an  explosion 
point  source.  However,  unlike  the  classical  explosion  model,  this  cyhndrical  cavity  source  does 
generate  SH-wave  energy  that  has  a  strong  azimuthal  amplitude  variation  (Figure  4).  Hence, 
even  though  the  S-wave  is  not  too  strong  at  the  source,  the  development  of  the  Lg  wavetrain  is 
still  efficient  enough  that  the  cavity  shape  has  a  detectable  influence  at  regional  distances.  We 
have  also  computed  3-D  finite  difference  synthetic  seismograms  for  an  earth  model  essentially 
the  same  as  the  one  described  in  Table  2,  except  with  a  randomly  irregular  Moho.  In  this  case, 
the  effect  of  the  cavity  shape  still  is  more  effective  in  generated  a  transverse  component  signal 
than  is  the  Moho  roughness. 

CONCLUSIONS  AND  RECOMMENDATIONS 

The  BEM  provides  a  useful  means  of  simulating  the  seismic  wavefields  generated  by  explosion 
sources  in  cylindrical  cavities.  Unlike  the  finite  difference  method,  which  requires  a  discretiza¬ 
tion  of  the  entire  earth  model,  an  implementation  of  the  BEM  requires  only  that  the  boundary 
and  the  cavity  be  divided  into  elements,  allowing  an  efficient  simulation  of  the  elastic  waves  at 
large  distances  from  the  source  cavity.  Results  for  cavities  in  a  hard  rock  situation  show  that 
though  there  are  significant  variations  in  amplitudes  of  P- waves  at  very  high  frequencies  (e.g., 
150  Hz),  at  frequencies  more  typical  of  regional  wave  propagation,  only  a  modest  change  in 
amplitude  is  predicted.  Likewise,  only  a  relatively  small  shear  wave  is  predicted.  An  explosion 
source  located  in  a  cylindrical  cavity  in  a  transversely  isotropic  formation  can  have  shear  wave 
amplitudes  even  lower  than  in  similar  isotropic  media.  Incorporation  of  these  explosion  sources 
into  discrete  wavenumber  models  of  regional  propagation  suggest  that  there  still  will  be  some 
azimuthal  variation  of  Lg  wave  amplitude  depending  on  propagation  direction  relative  to  the 
cavity  axis  (when  the  cavity  is  horizontal). 

Future  work  should  include  further  testing  of  these  theoretical  predictions  against  observed 
field  data,.  The  best  scenario  would  be  tests  against  data  collected  where  a  source  cavity 
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Figure  4:  Synthetic  seismograms  computed  for  the  earth  model  described  in  Table  2.  Receivers 
were  located  at  a  constant  offset  of  500  km  and  cover  the  full  range  of  azimuth  to  show  the 
effect  of  a  source  located  in  a  cylindrical  cavity  on  regional  wave  propagation.  Both  component 
plots  have  the  same  plotting  gain  applied  for  direct  comparison  of  amplitudes.  The  cavity  axis 
direction  corresponds  to  0°.  (A)  Vertical  component.  (B)  Transverse  component. 
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orientation  and  approximate  shape  is  known,  as  well  as  the  location  of  the  source  within  the 
cavity.  Another  valuable  task  is  further  consideration  of  simulations  where  we  include  equivalent 
moment  tensors  in  simulation  of  regional  wave  propagation  in  order  to  explore  the  effects  of  the 
explosion  sources  on  such  features  as  Lg  wave  amplitude.  The  overall  goal  is  to  better  develop 
our  understanding  of  near  source  effects  in  order  to  better  understand  the  influence  on  regional 
seismic  wave  propagation. 
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Parameter 

Background  Value  (GPa) 

Perturbed  Value  (GPa) 

C^ii 

35.2  GPa 

33.8  GPa 

Czz 

35.2 

26.0 

C'ls 

13.9 

10.3 

C44 

10.6 

9.53 

Cee 

10.6 

10.6 

P 

2.2  g/cm^ 

2.2 

Table  1:  Elastic  constants  for  the  background  medium  and  a  medium  with  aligned  cracks 
computed  using  the  theory  of  Hudson  (1980,  1981).  The  background  parameters  correspond 
to  P  and  S-wave  velocities  of  4000  m/s  and  2200  m/s,  respectively,  and  the  perturbed  values 

were  obtained  using  a  crack  density  of  0.05. 


Layer  thickness 

V, 

14 

Density 

1  km 

12 

22 

(half  space) 

5.3  km/sec 
6.1 

6.6 

8.1 

3.1  km/sec 
3.5 

3.8 

4.7 

2.4  g/cm^ 
2.6 

2.7 

2.8 

Table  2:  Crustal  velocity  model  used  to  compute  the  synthetic  seismograms  in  Figure  4. 
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